function out = getData(dataOption,startYear)

% The yield curve data: end of month data
loadYields     = load('GSW_yieldsMd_Mat3Md.mat');

% Remove the one-month and the three month interest rate
loadYields.yieldsMd_Mat3Md = [loadYields.yieldsMd_Mat3Md(:,1) loadYields.yieldsMd_Mat3Md(:,4:end)];

% Only include bond yields with max maturity of 120 months
select                     = loadYields.yieldsMd_Mat3Md(1,2:end) <= 120;
loadYields.yieldsMd_Mat3Md = loadYields.yieldsMd_Mat3Md(:,[1==1,select]);
% Change from discrete to continuous compounding - data is already in continous compounding
%loadYields.yieldsOut(2:end,2:end) =  log(1+loadYields.yieldsOut(2:end,2:end)/100);

% %CRSP risk-free rates for h = 1 and h = 3, in annualized terms
load('crspRet.mat');
crspRet1 = 12*log(1+crspRet(2:end,2));
crspRet3 = log(1+crspRet(2:end,3));
crspRet3 = 4*filter(ones(3,1),1,crspRet3);


% Loading all GSW yields (only needed for excess return regressions)
tmp = load('GSW_yieldsMd_Mat1Md.mat');
% TEST
%tmpMagnus = load('GSW_yieldsMd_Mat1Md_2016.mat');
%tmpMagnus.yieldsMd_Mat1Md(2:end,2:end) = log(1+tmpMagnus.yieldsMd_Mat1Md(2:end,2:end)/100);
%max(abs(tmpMagnus.yieldsMd_Mat1Md(2:end,120)-tmp.yieldsMd_Mat1Md(2:end-12,120)/100))


% Load macro variables
% tmpMacro   = xlsread('macroVariables_MacroFinanceModel.xlsx','Sheet1');
% macro.macroDates = x2mdate(tmpMacro(:,1));
% macro.LN   = tmpMacro(:,1+1);  %Ludvigson & Ng factor
% macro.CPI  = tmpMacro(:,1+2);  %Cieslak and Pavola factor for trend inflation
% macro.cpGAP= tmpMacro(:,1+3);  %Cooper & Priestley gap
% macro.PMI  = tmpMacro(:,1+4);  %PMI variable
% macro.NBER = tmpMacro(:,1+5);  %NBER recession indicator
% macro.ADS  = tmpMacro(:,1+6);  
% macro.CNAI = tmpMacro(:,1+7);  
% macro.RECPBROB = tmpMacro(:,1+8);  
% macro.coreCPI  = tmpMacro(:,1+9);
% macro.UGAP     = tmpMacro(:,1+10);
% macro.gdpCPI   = tmpMacro(:,1+11);    % inflation for the GDP deflator
% macro.potentOUT= tmpMacro(:,1+12);    % potential output as in Restud paper by Ang et al
% save('macroVariables_MacroFinanceModel.mat','macro')
macroLoad = load('macroVariables_MacroFinanceModel.mat');
macro     = macroLoad.macro;
%macro    = load('macroVariables_MacroFinanceModel_old.mat');

% Only include data from "startYear" and ending at 2016:12 (macro data does not go further)
% For yields
select1 = (str2num(datestr(loadYields.yieldsMd_Mat3Md(2:end,1),10)) >= startYear);
select2 = (str2num(datestr(loadYields.yieldsMd_Mat3Md(2:end,1),10)) <= 2016);
selectYield = select1.*select2 ==1;
loadYields.yieldsMd_Mat3Md = loadYields.yieldsMd_Mat3Md([1==1;selectYield],:);
loadYields.yieldsMd_Mat3Md(2:end,2:end) = loadYields.yieldsMd_Mat3Md(2:end,2:end)/100; %yields in base units

% For all yields
tmp.yieldsMd_Mat1Md = tmp.yieldsMd_Mat1Md([1==1;selectYield],:);
tmp.yieldsMd_Mat1Md(2:end,2:end) = tmp.yieldsMd_Mat1Md(2:end,2:end)/100; %yields in base units

% For macro (and CRSP yields)
select1 = (str2num(datestr(macro.macroDates(1:end,1),10)) >= startYear);
select2 = (str2num(datestr(macro.macroDates(1:end,1),10)) <= 2016);
selectMacro = select1.*select2 ==1;


%% Generating output
out.yieldsOut        = loadYields.yieldsMd_Mat3Md;
out.yields_h1        = tmp.yieldsMd_Mat1Md;
out.crspRiskFree_h1  = crspRet1(selectMacro,1);
out.crspRiskFree_h3  = crspRet3(selectMacro,1);
out.NBERrec          = macro.NBER(selectMacro,1);   

% T = sum(selectMacro);
% rho = 0.98;
% meanInfl = zeros(T,1);
% meanInfl(1,1) = macro.coreCPI(1,1);
% for t=2:T
%     meanInfl(t,1) = (1-rho)*macro.coreCPI(t,1) + rho*meanInfl(t-1,1);
% end
%plot([meanInfl/100,macro.CPI])


%out.infl             = (macro.coreCPI-0*meanInfl)/100;   %MMA = runHamiltonFilter(macro.coreCPI)
out.infl      = macro.coreCPI/100;
out.econ      = (macro.PMI(selectMacro,1)-44.5)/100;    % we use PMI: threshold value is 44.5

if dataOption == 1
    out.ygap             = -macro.UGAP/100;
elseif dataOption == 2
    % Using potential output to define the output gap
    out.ygap      = macro.potentOUT*100;
    %plot([macro.potentOUT*100,-macro.UGAP])
    %corr(macro.potentOUT*100,-macro.UGAP)
elseif dataOption == 3
    % Loading Cooper & Priestley output gap of the output gap
    out.ygap = macro.cpGAP(:,1);
    %plot([macro.cpGAP,-macro.UGAP])
    %corr(macro.cpGAP,-macro.UGAP)
elseif dataOption == 4
    % Loading the Ludvigson and Ng factor as the output gap
    out.ygap = macro.LN(:,1);
    %plot([macro.LN(:,1),-macro.UGAP])
    %corr(macro.LN(:,1),-macro.UGAP)
end 

end

